#!/usr/bin/env python

from scipy import comb

def basis_fun_eval(j,n,a,b,x):
    if(j > n or j < 0):
        return 0. # important for derivative evaluation
    else:
        return comb(n,j,exact=1)*(x-a)**j*(b-x)**(n-j)/(b-a)**n

def basis_fun_derivative(k,j,n,a,b,x):
	if (k > n): 
		print 'bernstein_basis_fun_derivative::Error: k > n!'
	elif (k == 1):
		return n/(b-a)*( basis_fun_eval(j-1,n-1,a,b,x) - basis_fun_eval(j,n-1,a,b,x) ) 
	else:
		return n/(b-a)*(basis_fun_derivative(k-1,j-1,n-1,a,b,x)-basis_fun_derivative(k-1,j,n-1,a,b,x))

def interpolant(beta,n,a,b,x):  
	suma = 0.
	for j in range(n+1):
		suma += beta[j]*basis_fun_eval(j,n,a,b,x)
	return suma
